Self-Contact and Instabilities in the Anisotropic Growth of Elastic Membranes 
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We investigate the morphology of thin discs and rings growing in the circumferential direction. 
Recent analytical results suggest that this growth produces symmetric excess cones (e-cones). We 
study the stability of such solutions considering self-contact and bending stress. We show that, 
contrary to what was assumed in previous analytical solutions, beyond a critical growth factor, no 
symmetric e-cone solution is energetically minimal any more. Instead, we obtain skewed e-cone 
solutions having lower energy, characterized by a skewness angle and repetitive spiral winding with 
increasing growth. These results are generalized to discs with varying thickness and rings with holes 
of different radii. 

PACS numbers: 05.45.-a, 46.70.Hg, 89.75.Da 
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Growing membranes are ubiquitous in living matter, 
producing a very rich variety of complex shapes. In an 
attempt to sustain external loads and to minimize inter- 
nal stresses, structures adopt new shapes. Non-uniform 
or restrained growth of surfaces inevitably leads to out of 
plane buckling. This can be evidenced in many examples 
such as the development of the folds of Mitochondria, the 
formation of plant leaves, the wrinkling of human skin or, 
in two dimensions, the packing of rod-like objects [lj. 

Flat discs growing under various conditions have re- 
ceived intense attention recently from experimental [2j |3] 
and theoretical sides 0HB]. Santangelo [4 presented 
a theoretical study to explain experiments on radially 
symmetric, isotropic growth of polymers with different 
swelling factors [2j[3]. In Dervaux et al. [5], potato-chip 
like shapes were found for anisotropic, circumferential 
growth using modified Foppl-von Karman equations [7]. 
This description is however only valid for small defor- 
mations and was extended in Ref. [6 to large deforma- 
tions in the limit of vanishing thickness. It was found 
that circumferential growth is equivalent to inserting a 
wedge into the disk, leading to so-called excess cones or 
" e-cones" . They are characterized by an excess angle 
ip e which is the angle of the wedge added to the disc. 
With increasing ip e , the symmetric e-cone will eventu- 
ally touch itself simultaneously with two perpendicular 
contact planes. From that situation on, the authors of 
Ref. [6 made the assumption that the e-cone remains 
symmetric, and derived the corresponding bending en- 
ergy. 

In this Letter, we show that this assumption is not 
correct when cp e is above a critical value. Instead, by 
taking surface self-contact into account, a family of new 
shapes emerges, namely skewed e-cone solutions with a 



varying contact plane angle. Our results are obtained 
from a sophisticated numerical thin shell model capa- 
ble of capturing large deformation and growth. We find 
the energetically most favorable shapes, discuss the de- 
pendence of stable morphologies on shell thickness, and 
highlight the analogy between disk and ring growth. 

The mechanics of thin discs is well described by the 
classical Kirchhoff-Love model. It follows by expansion 
of the full 3d elastic problem in the thickness ft, and by 
assuming negligible transverse shear through the thick- 
ness [8-10 . The elastic energy is then entirely given by 
the deformation of the disc's middle-surface Q, and con- 
tains a stretching and a bending part: 
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Here, Y is the elastic modulus, v Poisson's ratio and 
d£l the infinitesimal area element. The stretching and 
bending energy densities U s and [/& are given by 

U 8 = (eii+e 22 ) 2 -2(l-^)(e 11 e 22 -e? 2 ) and (2) 
U b = («n + ft 22 ) 2 -2(1 - v)(k,uk,22 - ^12) • ( 3 ) 

€ij denote the in-plane strains and Kij the components of 
the curvature tensor. The first term in Ub is the squared 
mean curvature, while the second term accounts for 
Gaussian curvature. Without growth, e = 1/2(F T F — 1), 
where F = Vx is the deformation gradient of the map 
X that transforms from the reference to the deformed 
configuration, see e.g. Ref. [TT] . Growth is added by 
decomposing F multiplicatively into F = AG jT2j [13]. 
G is the growth tensor describing the change of mass 
and A is the elastic tensor that ensures compatibility (no 
overlap) and continuity (no cavitation) of the body. We 
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FIG. 1: Numerical equilibrium solution for the touching prob- 
lem of e-cones: Increasing the excess angle cp e of a symmetric 
e-cone (a), the four wings start to touch, initially with angle 
7 = 7r/2 between the two touching planes (b). For larger ip e , 
the angle changes (c) and a skewed e-cone emerges (d), where 
the contacting wings have slided against each other. A paper 
model is shown in (e) for comparison. Colors in (a) and (d) 
represent the bending energy density Ub> Red color in (b) 
and (c) denotes points with self-contact. Simulation parame- 
ters are Y = 10 4 Pa, v = 0.3, R = lm, h = 0.02 m (lenghts 
are rescaled by R in the main text). Values of cp e for (a)-(e) 
are (3.73, 7.45, 7.6, 47T, 4tt). Dashed lines in (a) denote points 
P c of vanishing curvature. For ip e = 3.73, their tangents are 
perpendicular to the gray-shaded plane. 



further assume: (i) Growth is slow compared to the elas- 
tic timescale, which is typically true in nature [13]. Thus, 
the disc is always in elastic equilibrium, (ii) There ex- 
ists a stress-free reference state, i. e. the flat disc with 
G — 1. (Hi) The elastic response depends only on A. 
The in-plane strains then become e = 1/2(A T A — 1) and 
k is modified similarly. In the present model, we only 
consider surface expansion, i. e., the thickness remains 
constant. In terms of the excess angle <^ e , G in polar 
coordinates then becomes 



G(r,<p) 
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We search for configurations that minimize Eq. (fil) by 
discretizing the system using the finite element method. 
The curvature integral in Eq. (fil requires at least con- 
tinuity of first derivatives. To this end, we employ a 
recently developed method based on C 1 -continuous sub- 
division finite elements (SDFEs) [14]. From the discrete 
Lagrangian formulation of Eq. (fil the elastic forces f^ 
are derived for all mesh nodes z, and the equations of 
motion are integrated in time. For numerical stability 
and equilibration, we add an external viscous damping 



term of the form rjVi to the equations of motion, with v^ 
the velocity of node i. r] was chosen subcritical for all vi- 
bration modes and did not have influence on the results. 
To validate our model, we use the non-linear ABAQUS 
software package with implicit solver and 4 node gen- 
eral shell elements (S4) including transverse shear [15]. 
Contact is modeled in ABAQUS via surface- surface con- 
tact and a linear force law while in our implementation, 
point-triangle contacts are subject to a quadratic repul- 
sion force. Note that in tangential direction contacts are 
frictionless unless explicitly mentioned. 

The simulation starts from a flat disc with constant 
thickness h, with h between 0.001 and 0.04. The excess 
angle ip e is increased in small steps by expanding all ele- 
ments in the tangential direction of a centered cylindrical 
material coordinate system. Assuming slow growth, the 
equilibrium solution is obtained after every such step (in 
SDFEs, we wait until the kinetic energy becomes negligi- 
ble). In SDFE simulations a small vertical random noise 
is initially imposed on the displacements to break the 
planar symmetry of the flat disk. In ABAQUS, this is 
automatically provided by numerical noise. In dynamic 
simulations via ABAQUS, we first observe high buckling 
modes of the e-cone type, and for h small enough even 
superimposed radial buckles. Due to the bending rigid- 
ity these soon damp out into the lowest energy mode 
with two folds (n=2), having a shape similar to a potato 
chip as shown in Fig. [I] (a). The growth continues sta- 
ble in the n=2 mode. The analytical solution predicts 
contact between the two folds at (p^ lss = 7.08 [6 . We 
find excellent numerical agreement, with the first con- 
tact occurring between cp e = 7.12 ± 0.01 for h = 0.04, 
and <p e = 7.074 ± 0.005 for h = 0.001. 

The analytical solution suggests an increase in bend- 
ing energy up to the theoretical value of the solution 
with three wings (n=3), followed by a transition into 
this mode [6]. We obtain a flattening of the contact 
zones maintaining symmetry in two planes (Fig. [11(b)). 
Shortly after reaching (/?Jf ss , however, the contact planes 
start to counter rotate (Fig. fl] (c)), leading to a reduc- 
tion of curvature of one wing at the expense of the other. 
Consequently a skewed e-cone solution with a flat disc- 
like part and two wings forms as shown in Fig. [I] (d). 
This transition is obtained for all shell thicknesses con- 
sidered. If growth continues, we simply keep on adding 
additional windings to both sides of the flattened region. 
In the following we address the nature of the transition, 
the composition of elastic energies, and their dependence 
on system parameters. 

The analytical results in Ref. [6 are based on discs 
with vanishing thickness. As a consequence, the e-cone 
shape has zero Gaussian curvature everywhere, except at 
the point-like apex. By introducing a finite thickness, it 
is intuitively clear that the singularity at the apex will be 
flattened in order to reduce bending energy and because 
the maximal curvature is bounded by the thickness. Con- 
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FIG. 2: Normalized bending energy Eb in dependence on the 
disc thickness h for the symmetric e-cone, compared to the 
analytical solution with h = [6]. Inset: The position of the 
maxima of Eb, denoted (£™ ax , approaches the analytical value 
ipt according to (^ ax - ip* e ) oc h a with a = 0.45 ± 0.03. 



sequently, the Gaussian curvature cannot vanish close to 
the apex leading to local stretching. We call this region 
the core of the e-cone and denote its radius by R c . We 
numerically calculated the e-cone equilibrium solutions 
for different thicknesses and excess angles and found that 
R c oc h, a behavior similar to the core size scaling of de- 
velopable cones reported in [HJ [TBI HZ] - In Fig. [2] we show 
the rescaled bending energy E b = j Q U b /A dVt, where 
A is the total area of the e-cone [18] . The figure shows 
that the bending energy is indeed strongly reduced for 
thick shells. To quantify the convergence towards the 
analytical solution, we consider the excess angle (^™ ax 
where E b reaches its maximum. We find that (^™ ax con- 
verges to the analytical maximum tp* = 2.57 according 
to (^ ax - ip* e ) (x h a with a = 0.45 ± 0.03 (Fig. [2| inset). 

For the n = 2 mode at (/?JP SS , the wings touch and even- 
tually form a skewed e-cone. In order to analyze the sta- 
bility of the symmetric solution, we performed two types 
of simulations: First, we "stabilized" the symmetric so- 
lution using static friction between the areas in contact 
and increased cp e to 10. This way, we can fix the e-cone 
solution, which would be unstable without friction. At 
cp e = 10, the contact friction was removed and the sys- 
tem relaxed towards the skewed e-cone. We then followed 
the skewed solution by decreasing the excess angle from 
cp e = 10 back to cp^ lss . During this procedure, we mea- 
sured the bending energy, which is the relevant part of 
the total energy, shown in Fig. [3] We note that in the 
range of ip e investigated, the next higher n = 3 mode 
(dotted line) is energetically above the symmetric e-cone 
(bold line), despite the fact that it is a configuration that 
does not involve any contact (for n = 3, the wings would 
touch only for (p e = 13.3). Fig. [^furthermore shows that 
above a critical (/?g nt , the skewed solution provides the 
lowest energy, whereas the symmetric e-cone is favored 



FIG. 3: Bending energies for the n — 2 contact problem with 
h — 0.02: The symmetric n — 2 e-cone (bold line) is stable 
until ip e = 8, when the solution changes into the skewed e- 
cone (thin line). For comparison, the unphysical situation 
with self- intersections is also shown (dashed line), as well as 
the n — 3 solution (dotted line). The inset illustrates the 
influence of the thickness in the transition region. Bold lines 
are again the symmetric e-cones, and thin lines correspond 
to the skewed equilibrium solution. The critical excess angles 
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where e-cones loose stability are marked with black dots. 



in the regime (p^ lss < cp e < (fl rit . The influence of the 
finite disc thickness on the energy and stability is illus- 
trated in the inset of Fig. [3] Remarkably, even though 
(/?g lss is practically identical for all thicknesses considered, 
we find a significant dependence of (/?g nt on h (black dots) . 
Finding the minimum of the nonlinear problem posed 
by Eq. (IT]) with self-contact is nontrivial and numerically 
delicate. We therefore compared our numerical findings 
to experiments by measuring the angle 7 between the 
touching planes (cf. Fig. IT]) for h=0.005. In simulations, 
this was done by averaging over the contact forces in each 
of the two contact zones, resulting in an approximation 
for the contact plane normals from which the angle 7 is 
obtained. Experimentally, we used two layers of paper 
glued together to obtain a homogenous bending stiffness. 
The angle was obtained from digital camera images taken 
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FIG. 4: Left: Dependence of the touching plane angle 7 
(cf. Fig. [l]) on (f ej determined numerically (circles) and ex- 
perimentally (squares) for h = 0.005. Right: Bending energy 
ratio / = Eb&i ouch )/E h (Ve = 3.73) • (27r + ^ e ouch )/(27r + 3.73) 
between skewed e-cones and the symmetric one at cp e — 3.73. 
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FIG. 5: Numerical equilibrium solution (left) for a growing 
ring compared to a paper model (right) at ip e = 4tt. Colors 
represent the bending energy density [/&■ Simulation param- 
eters are identical to Fig.[l] with Ri = 0.5. 



^e 



Ri = 0.1 Ri = 0.25 Ri = 0.5 



7.3 ±0.05 7.8 ±0.05 8.45 ± 0.05 



TABLE I: Dependence of (p e lss on the inner radius Ri of a 
ring for h = 0.02. 



that were not predicted by previous analytical solutions, 
emerge and stabilize for sufficient growth. By energeti- 
cal arguments based on symmetric deformation modes, 
we showed, that the asymmetric mode has lower bending 
energy, and is therefore preferred. Moreover, we found 
that tangent ially growing rings of various widths behave 
analogous to discs. Good agreement with experiments 
suggests that the skewed solutions are energetically op- 
timal, but a rigorous proof remains an open question. 
The asymmetries arising from circumferential growth dis- 
closed in the present work constitutes a novel step in un- 
derstanding the shapes of membranes and sheets. One 
can generalize our work to non-linear material behavior 
and stress or strain dependent growth laws in order to 
tune to more specific applications in Nature. 

This work was supported by Grant No. TH-06 07-3 
of the Swiss Federal Institute of Technology Zurich and 
FUNCAP. The authors would like to thank J. Guven for 
helpful discussions. 



from the top view of the e-cone. The results in Fig. [4] 
(left) show good agreement of simulations and experi- 
ments, confirming the validity of the model fc\l and the 
numerical equilibrium solution. 

Moreover, for values of ip e > 10 the skewed e-cone so- 
lution approaches a shape which can be approximated 
by a flat disk with two connecting loops (see Fig. Hid). 
Both loops put together are described remarkably well 
by the symmetric free e-cone shown in Fig. [I] a: From 
the shape equations of the symmetric e-cone it can be 
derived that for cp e = 3.73, the tangent in the points P c 
of vanishing curvature is perpendicular to the plane in 
which the P c lie. These points are in fact equivalent to 
the points of contact of the skewed e-cone. The described 
approximation yields an estimate for the bending energy. 
One obtains for a skewed e-cone with excess angle ^ouch. 
E b ((pi™ ch ) » E b (<p e = 3.73) • (2tt + 3.73)/(2tt + (^ ouch ). 
In Fig. p] (right), we plot the ratio of the two energies 
Eb(<Pe)/M<Pe = 3.73) • (2tt + ^ e ouch )/(27r + 3.73), 
which is indeed close to the predicted value (solid line). 
Understandably, the quality of the approximation de- 
creases with larger h, as the change in bending energy 
of the e-cone core region is not taken into account. 

We finally turn to the study of rings with inner radius 
Ri. The same simulation procedure was applied with 
fixed thickness h = 0.02. We found that rings, like discs, 
loose stability after contact occurs and attain the skewed 
solution shown in Fig. ^1 Interestingly, (^jj lss strongly 
depends on Ri , as shown in Table [TJ In terms of energy 
and stability, rings behave similar to discs. 

We presented a realistic simulation of circumferential 
growth of discs and rings with finite thickness. By in- 
cluding self-contact, asymmetric, skewed morphologies, 
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